########################
###19.fe_single_int.R#############
########################

pols <- readRDS('pols_cwg')

#####Fixed Effects Models#####

library(plm)
library(pcse)
library(lmtest)
library(dplyr)

#####gcov:gini#####

gg1 <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc1 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc1 + 
             scale_galtansdw +
             scale_enep +
             scale_unemp,
           data = pols, 
           na.action = na.omit, 
           model = "within")

coeftest(gg1, vcov = function(x) vcovBK(x, cluster="group"))

gg2 <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc2 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc2 + 
             scale_galtansdw +
             scale_enep+
             scale_unemp,
           data = pols, 
           na.action = na.omit, 
           model = "within")

coeftest(gg2, vcov = function(x) vcovBK(x, cluster="group"))

gg3 <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc3 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc3 + 
             scale_galtansdw +
             scale_enep+
             scale_unemp,
           data = pols, 
           na.action = na.omit, 
           model = "within")

coeftest(gg3, vcov = function(x) vcovBK(x, cluster="group"))

gg4 <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc4 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc4 + 
             scale_galtansdw +
             scale_enep+
             scale_unemp,
           data = pols, 
           na.action = na.omit, 
           model = "within")

coeftest(gg4, vcov = function(x) vcovBK(x, cluster="group"))

gg5 <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc5 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc5 + 
             scale_galtansdw +
             scale_enep+
             scale_unemp,
           data = pols, 
           na.action = na.omit, 
           model = "within")

coeftest(gg5, vcov = function(x) vcovBK(x, cluster="group"))

##### gini*ggini Point Estimates#####

fu1 <- as.data.frame(summary(gg1, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) 

fu2 <- as.data.frame(summary(gg2, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(gg3, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(gg4, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(gg5, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Unemployment", "Gini*Income Dif.")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini", "Gini*Income Dif.", "GALTAN Polarization","ENEP", "Unemployment"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))

library(ggplot2)

dot.table.fu %>% ggplot(aes(x = Predictor, y = pe)) + 
  geom_point() + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed') +
  geom_linerange(aes(ymin = lwr, ymax = up))+
  coord_flip() +
  labs(title = "Point Estimates of Model Regressing Predictors on Economic Polarization (n = 83)", y = "Point Estimates (95% Confidence Interval)", x = "Predictor")+
  theme_bw()

rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))


#Model Fit Statistics#

fit1i <- c(summary(gg1)$r.squared, summary(gg1)$fstatistic$statistic, summary(gg1)$fstatistic$p.value)
fit2i <- c(summary(gg2)$r.squared, summary(gg2)$fstatistic$statistic, summary(gg2)$fstatistic$p.value)
fit3i <- c(summary(gg3)$r.squared, summary(gg3)$fstatistic$statistic, summary(gg3)$fstatistic$p.value)
fit4i <- c(summary(gg4)$r.squared, summary(gg4)$fstatistic$statistic, summary(gg4)$fstatistic$p.value)
fit5i <- c(summary(gg5)$r.squared, summary(gg5)$fstatistic$statistic, summary(gg5)$fstatistic$p.value)


fit.statsi <- rbind(fit1i, fit2i, fit3i, fit4i, fit5i)

print(mean(fit.statsi[,1])) #R-Squared
print(mean(fit.statsi[,2])) #Adjusted R-Squred
print(mean(fit.statsi[,3])) #F-Statistic
print(mean(fit.statsi[,4])) #p-value of F-Statistic

####Marginal Effect Plot####

library(MASS)
library(gridExtra)
#####Figure 2: Sorting on X-axis, salience held constant in each panel#####

pols_rug <- readRDS('pols_bjps') %>%
  dplyr::select(scale_econsdw,scale_dispgini, 
                scale_ggini_inc1,
                scale_socioeconomicsal,
                scale_galtansdw ,
                scale_enep ,
                scale_unemp) %>%
  na.omit
#simulate 1000 draws from a multivariate random normal distribution based on the model#
test <- as.data.frame(mvrnorm(1000, mu =coef(gg1), Sigma =vcovHC(gg1)))

#vector of sorting from observed min to observed max#
nums <- seq(from = min(pols$scale_ggini_inc1, na.rm = T), to = max(pols$scale_ggini_inc1, na.rm = T), length.out = 1000)

#create data frame for storing output#
marg <- data.frame(x = nums)

#Simulate marginal effect based on draws from the model#
for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1]
  }
}

#Vectors to be filled#
sort <- NA
pe <- NA
lb <- NA
ub <- NA


#Make CI based on empirical quantiles

for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

#Combine into a singe dataset

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

#create plot. 

sort_int <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-2, 3.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "", x = "Partisan Sorting by Income", y = "Marginal Effect of Disposable Income Inequality on Party Polarization") +
  theme_bw()+
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

#####sal:gini#####

sal1 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc1 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_socioeconomicsal + 
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(sal1, vcov = function(x) vcovBK(x, cluster="group"))

sal2 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc2 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_socioeconomicsal + 
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(sal2, vcov = function(x) vcovBK(x, cluster="group"))

sal3 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc3 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_socioeconomicsal + 
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(sal3, vcov = function(x) vcovBK(x, cluster="group"))

sal4 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc4 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_socioeconomicsal + 
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(sal4, vcov = function(x) vcovBK(x, cluster="group"))

sal5 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc5 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_socioeconomicsal + 
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(sal5, vcov = function(x) vcovBK(x, cluster="group"))

##### gini*sal Point Estimates#####

fu1 <- as.data.frame(summary(sal1, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) 

fu2 <- as.data.frame(summary(sal2, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(sal3, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(sal4, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(sal5, vcov = function(x) vcovBK(x, cluster="group"))$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Unemployment", "Gini*Salience")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini", "Gini*Salience", "GALTAN Polarization","ENEP", "Unemployment"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))

library(ggplot2)

dot.table.fu %>% ggplot(aes(x = Predictor, y = pe)) + 
  geom_point() + 
  geom_hline(aes(yintercept = 0), linetype = 'dashed') +
  geom_linerange(aes(ymin = lwr, ymax = up))+
  coord_flip() +
  labs(title = "Point Estimates of Model Regressing Predictors on Economic Polarization (n = 83)", y = "Point Estimates (95% Confidence Interval)", x = "Predictor")+
  theme_bw()

rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))

#Model Fit Statistics#

fit1h <- c(summary(sal1)$r.squared, summary(sal1)$fstatistic$statistic, summary(sal1)$fstatistic$p.value)
fit2h <- c(summary(sal2)$r.squared, summary(sal2)$fstatistic$statistic, summary(sal2)$fstatistic$p.value)
fit3h <- c(summary(sal3)$r.squared, summary(sal3)$fstatistic$statistic, summary(sal3)$fstatistic$p.value)
fit4h <- c(summary(sal4)$r.squared, summary(sal4)$fstatistic$statistic, summary(sal4)$fstatistic$p.value)
fit5h <- c(summary(sal5)$r.squared, summary(sal5)$fstatistic$statistic, summary(sal5)$fstatistic$p.value)


fit.statsh <- rbind(fit1h, fit2h, fit3h, fit4h, fit5h)

print(mean(fit.statsh[,1])) #R-Squared
print(mean(fit.statsh[,2])) #Adjusted R-Squred
print(mean(fit.statsh[,3])) #F-Statistic
print(mean(fit.statsh[,4])) #p-value of F-Statistic

#####Make Marginal Effects Plot


#simulate 1000 draws from a multivariate random normal distribution based on the model#
test <- as.data.frame(mvrnorm(1000, mu =coef(sal1), Sigma =vcovHC(sal1)))

#vector of sorting from observed min to observed max#
nums <- seq(from = min(pols$scale_socioeconomicsal, na.rm = T), to = max(pols$scale_socioeconomicsal, na.rm = T), length.out = 1000)

#create data frame for storing output#
marg <- data.frame(x = nums)

#Simulate marginal effect based on draws from the model#
for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] 
  }
}

#Vectors to be filled#
sort <- NA
pe <- NA
lb <- NA
ub <- NA


#Make CI based on empirical quantiles

for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

#Combine into a singe dataset

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

#create plot. 

sal_int <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "", x = "Salience of Economic Issues", y = "Marginal Effect of Disposable Income Inequality on Party Polarization") +
  theme_bw()+
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

library(gridExtra)

fe_single_interactions <- grid.arrange(sort_int, sal_int, nrow = 1, top = "Marginal Effect of Income Inequality from Fixed Effects Models", left = "", bottom = "")

#ggsave("figure_A9.pdf", plot = fe_single_interactions, width = 11, height = 5, units = 'in')
